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ABSTRACT 

A molecule traveling in a realistic propagation environment 
can experience stochastic interactions with other molecules 
and the environment boundary. The statistical behavior of 
some isolated phenomena, such as dilute unbounded molec¬ 
ular diffusion, are well understood. However, the coupling of 
multiple interactions can impede closed-form analysis, such 
that simulations are required to determine the statistics. 
This paper compares the statistics of molecular reaction- 
diffusion simulation models from the perspective of molecu¬ 
lar communication systems. Microscopic methods track the 
location and state of every molecule, whereas mesoscopic 
methods partition the environment into virtual containers 
that hold molecules. The properties of each model are de¬ 
scribed and compared with a hybrid of both models. Simula¬ 
tion results also assess the accuracy of Poisson and Gaussian 
approximations of the underlying Binomial statistics. 

1. INTRODUCTION 

The prevalence of using molecules to communicate in bi¬ 
ological systems (see Ch. 16]) has recently attracted the 
attention of the research community to adapt the princi¬ 
ples of molecular communication (MC) for new applications 
to transmit arbitrary amounts of information in environ¬ 
ments where conventional methods of communication might 
be hazardous or impractical; see [TT]. One MC method, free 
diffusion, is attractive because it does not require additional 
infrastructure in the propagation medium. Free diffusion is 
effectively a random process where a molecule collides with 
other molecules in a fluid environment. 

The behavior of any one molecule in a realistic propaga¬ 
tion environment is unlikely to be characterized by diffu¬ 
sion alone. Other potential phenomena include bulk fluid 
flow, collisions with the environment boundary, and chemi¬ 
cal reactions either throughout the environment or in a lo- 


Permission to make digital or hard copies of all or part of this work for personal or 
classroom use is granted without fee provided that copies are not made or distributed 
for profit or commercial advantage and that copies bear this notice and the full citation 
on the first page. Copyrights for components of this work owned by others than the 
author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or 
republish, to post on servers or to redistribute to lists, requires prior specific permission 
and/or a fee. Request permissions fi'om Permissions@acm.org. 

NANOCOM '15, September 21 - 22 2015, Boston, MA, USA 

Copyright is held by the owner/author(s). Publication rights licensed to ACM. 

ACM X-XXXXX-XX-X/XX/XX. 


cal region. Generally, these phenomena contribute to the 
stochastic behavior of any single molecule. 

We noted in 13 that communications analysis requires 
the form and the statistics of the end-to-end channel im¬ 
pulse response, i.e., the time-varying signal observed at the 
receiver given that molecules are released at some instant by 
the transmitter. The response can then be used to derive 
the received signal for any modulation scheme. Analytical 
models for some isolated processes are known, such as for 
molecular diffusion; see [^. However, when multiple inter¬ 
actions are present, their impact is coupled and this can 
impede closed-form theoretical analysis. Often, simplifying 
assumptions are made and specific geometries are studied to 
facilitate analysis. For example, we analyzed an unbounded 
environment with diffusion, bulk fluid flow, and molecule 
degradation in . A closed-form time domain channel im¬ 
pulse response was derived, but this was in the absence of 
any local chemical reactions (such as at the receiver). Gen¬ 
erally, we may need to rely on numerical methods or simu¬ 
lations to determine the channel statistics. 

Simulation methods for molecular behavior can range in 
scale from molecular dynamics models (such as that used 
in LAMMPS |14| ), which account for all interactions be¬ 
tween all individual molecules (including solvent molecules 
in a fluid), to continuum models (such as that used in COM- 
SOL Multiphysics [^) where no individual molecules are de¬ 
scribed. Two common “intermediate” models that tend to 
be suitable for the study of reaction-diffusion environments 
are microscopic and mesoscopic models. Both of these mod¬ 
els treat the solvent in a fluid as a continuum and focus on 
the behavior of solute molecules. 

Microscopic simulators such as the Smoldyn simulator 
track the coordinates and behavior of each solute molecule; 
see [^. Mesoscopic simulators partition the system into vir¬ 
tual containers and track the number and type of solute 
molecules in each container. If molecular concentrations in 
each container are homogeneous, then a mesoscopic simu¬ 
lation can accurately capture the behavior of the system; 
see . However, the assumption of homogeneity can place 
severe constraints on the size of virtual containers; see 
|15| . A microscopic model has better spatial accuracy, but 
the advantages of a mesoscopic model include easier imple¬ 
mentation of complex chemical reactions and better compu¬ 
tational efficiency as the system dimensions grow. 

From a MC perspective, we are ultimately interested in 



the accuracy of the statistics at the receiver. We may need 
to simulate the system many thousands of times to com¬ 
pile the receiver statistics. The behavior of the total system 
is not as important, so we are motivated to improve com¬ 
putational efficiency in regions that are not critical to the 
receiver statistics, i.e., are sufficiently far from the commu¬ 
nication link. In [^, we combined two schemes towards 
this goal. In the first scheme, also described in , the local 
accuracy is adjusted by using mesoscopic containers of differ¬ 
ent sizes; generally, larger subvolumes are less accurate but 
more computationally efficient. In the second scheme, the 
environment is partitioned into microscopic and mesoscopic 
regimes, thus providing additional flexibility in the trade¬ 
off between local accuracy and computational complexity. 
These hybrid schemes have been proposed in papers includ¬ 
ing [T^[8l||. 

In this paper, we use simulations to study the accuracy 
of both the average time-varying channel response and the 
statistics of the response at specific times. Unlike in [13| , 
where we observed the channel response due to diffusion 
only, here we study the channel response and statistics of dif¬ 
fusion, a first-order chemical reaction, and a simple reaction- 
diffusion scenario. We gain insight into how aggressively we 
can reduce the computational complexity of the simulation 
environment without compromising the accuracy at the re¬ 
ceiver. A formal study of computational complexity is left 
for future work, but preliminary results were shown in [13| . 

We also compare the suitability of the Poisson and Gaus¬ 
sian approximations of the Binomial distribution when rep¬ 
resenting the cumulative distribution function (CDF) of re¬ 
ceiver observations made at a specihc time. These approxi¬ 
mations are commonly used in communications analysis be¬ 
cause they are less computationally intensive, and were also 


recently assessed for MC systems in 17 


The rest of this paper is organized as follows. The un¬ 
derlying physical model and channel statistics are described 
in Section The microscopic, mesoscopic, and hybrid sim¬ 
ulation models are briefly defined in Section Simulation 
results to compare the models are presented in Section 
Section concludes the paper. 


2. PHYSICAL MODEL 

In this section, we describe the common physical model 
that the simulation models represent. We present the ex¬ 
pected channel responses for the scenarios that we will sim¬ 
ulate, and discuss the statistics of those responses. 

The environment is a bounded two-dimensional fluid “vol¬ 
ume” V with a reflective boundary; generally, it could be 
absorbing or reactive if there are local chemical reactions at 
the boundary. There is a single molecular species, labeled 
molecule A, with constant diffusion coefficient D. Using a 
constant D implies that the A molecules are dilute. 

We consider two diffusive scenarios where the expected 
channel response at the receiver (RX) and the corresponding 
statistics are (at least approximately) known so that we can 
focus on assessing the accuracy of the simulation models. 
In the first scenario, we distribute molecules uniformly over 
V and observe the number of molecules present in Vrx, a 
subset of V. If A molecules are distributed, then the 
number of molecules expected in Vrx, URxit), is constant 
and equal to 

( 1 ) 


In the second scenario, we release (i.e., transmit) N A 
molecules in a small area near the center of V and observe 
(i.e., receive) the number of molecules present in another 
small area also near the center of V. If V is large enough to 
model as infinite, and the transmitter (TX) and RX regions 
are small enough to model as points, then from [Zl Eq. (3.4)] 
we can write Urx (t) as 


URx{t) = 


i-RX 


AnDt 


exp 



( 2 ) 


where fiRx is the area of the RX, and d is the distance be¬ 
tween the centers of the TX and the RX. Eq. § is accurate 
if hRx ^ d and d^ <C V. 

An A molecule can also degrade according to the first- 
order chemical reaction A —>■ 0, where k is the reaction rate 
constant in s“^. If there are N A molecules in the system 
at time t = Os, and the RX is all of V, then the number of 
molecules expected to remain at time t > 0, is § Eq. (9.7)] 

URx(t) = A^exp (-fct), (3) 


such that each molecule has a probability of exp {—kt) of 
remaining at time t, and we can account for degradation in 
the diffusive scenarios by scaling Q or (§ by exp (—kt). 

Now consider the statistics of the channel responses at a 
specihc instant t. Assuming no knowledge of one molecule’s 
location or whether it has been degraded after time t > 
0, then whether that molecule is in the RX at time t is 
the outcome of an independent trial; see [16[ Ch. 5.1]. For 
general N, there is one trial for each molecule. The number 
of molecules observed in the RX at time t, URx{t), is the 
number of “successful” trials. Thus, URx{t) for some t is a 
Binomial random variable, where the probability of success 
of each trial is equal to URx{t) with N = 1. 

By knowing URx{t), we can compare the empirical CDF 
of each simulation model with the Binomial CDF. We can 
also assess the Poisson and Gaussian approximations of the 
Binomial CDF. From |16| Ch. 5.2], the Poisson approxima¬ 
tion should be accurate when N is “large” and URx(t) for 

= 1 is “small”. From [16[ Ch. 5.5] and the Central Limit 
Theorem, the Gaussian approximation should be accurate 
for sufficiently large N. 


3. SIMULATION MODELS 

In this section, we summarize the simulation models that 
we will assess in Section [d] 

3.1 Microscopic Model 

In the microscopic model, the environment U is a single 
container Vm and there is a constant time step AIm. For 
each time step, the coordinates of every A molecule are up¬ 
dated by adding a random displacement ny/2DAtM to each 
dimension, where n is an independent normal random value 
with mean 0 and variance 1. Any molecule that ends up out¬ 
side of Um is reflected against the boundary of Vm. A given 
molecule is degraded during the time step and removed from 
17 if u > exp(—/ cAIm), where u is an independent uniform 
random value between 0 and 1. 

3.2 Mesoscopic Model 

In the mesoscopic model, V is partitioned into virtual sub¬ 
volumes or containers. We track the number of A molecules 
within each subvolume and assume that each subvolume is 


URx{t) = NVrxIV. 













homogeneous. Mesoscopic simulations are described as a 
series of “events”. A diffusion event is the transition of a 
molecule between adjacent subvolumes, and a degradation 
event is a decrement of the number of molecules in one sub¬ 
volume. Every event is assigned a propensity, a, which de¬ 
termines the probability of that event occurring next. In 
this paper, we summarize how the relevant propensities are 
calculated. Due to space, the reader is referred to or to 
works such as for details on how to use the propensities 
to simulate an event sequence. 

Diffusion propensities describe the expected transition 
rate of molecules between adjacent subvolumes. Following 
[13| , we consider square subvolumes that could have different 
sizes in order to adjust the local computational complexity. 
The transition rate from a square subvolume of width hi to 
one of width hj, where the overlap of their adjacent faces is 
length ho < mm{hi, hj}, is found to be 13 Eq. (9)] 
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where Ui is the number of molecules currently in subvolume 
i. The propensity of a chemical reaction describes the ex¬ 
pected frequency of the reaction. For our first-order degra¬ 
dation, the propensity Orxn.i of the reaction in subvolume i 
is given by Eq. ( 6 )] as flrxn.i = kUi. 


3.3 Hybrid Model 

In the hybrid model, V is partitioned into a microscopic 
regime Vm and a mesoscopic regime Vs. Both regimes are 
treated independently, as previously described, until there 
are molecules that transition from one regime to the other. 
For simplicity, as i n [13] , we adopt the simplified transition 
rules described in [ 10 ] : 

Vs to Vm'- a source subvolume in Vs must be along the 
boundary with Vm, and it has a mirror “imaginary” subvol¬ 
ume of the same size in Vm- A molecule leaving Vs to enter 
Vm is placed at random within the mirror subvolume and 
then treated as an individual molecule in Vm- 

Vm to Vs^ When a molecule in Vm is identified to have 
entered Vs, then we add that molecule to the subvolume 
along the boundary with Vm that is closest to the molecule’s 
new location. 


Figure 1: Simulation environment V drawn to scale. 
The inner region Vi has width 30 /rm and height 
15 pm. The TX and RX are squares of width 3 pm, 
are placed in the middle of Vi, and are separated 
by a distance of 15pm (center to center). The mid¬ 
dle region V 2 (in grey) surrounds Vi and has width 
90 pm. The outer region V 3 (in black) surrounds V 2 , 
has width 180 pm, and has a reflective outer bound¬ 
ary. The partitioning of the regions of V shown here 
is an example of the HYB model in Table [l] 


Table 1: System partitioning models, unless oth¬ 
erwise specified. When a region is mesoscopic, all 
subvolumes in that region have the specified width. 
V 2 for a HYB model is mesoscopic if /i 2 is defined. 
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4. NUMERICAL RESULTS 

In this section, we present simulation results to assess the 
accuracy of the simulation models to generate the channel 
response and the channel statistics. We consider the envi¬ 
ronment defined in Fig. (or a subset of it) for all of our 
simulations. The coefficient of diffusion is D = 10~®^. 
The environment V is partitioned into 3 regions. Each re¬ 
gion is modeled as either microscopic or mesoscopic. Unless 
otherwise specified, the model labels used in the simulation 
figures are described in Table Namely, we consider micro¬ 
scopic (MICRO), mesoscopic (MESO), multi-scale MESO 
(MESO-MS) and hybrid (HYB) partitioning models, where 
both the MESO-MS and HYB models are less accurate in 
V 2 and/or V 3 because the communication link is in Vi. 

We describe a series of 5 simulations as summarized in Ta¬ 
ble [I The first two simulations are uniform diffusion tests. 
For the first test, the system consists of region Vi only (this 
is the only test where molecule motion is restricted to Vi). 
1000 molecules are initialized over all of Vi and we observe 


the number of molecules present in one half of the region 
after t = 5s, i.e., 500 molecules are expected. The empirical 
CDF for each partitioning model, compiled over 10^ real¬ 
izations, is presented in Fig. The empirical CDFs of the 
simulation models all match the Binomial CDF, including 
the MESO-MS model where here we observe the number of 
molecules in a subvolume with width 15 pm and the rest of 
region is partitioned into subvolumes with width 1 pm. The 
Gaussian approximation is very close to the Binomial CDF 
whereas the Poisson approximation is not, since the under¬ 
lying trial success probability (i.e., the probability that a 
given molecule is observed) is 0.5, which is very high. 

In the second simulation, we perform a uniform diffusion 
test where we initialize 10800 molecules over all of V and we 
observe the number of molecules present in Vi after t = 20 s. 
From 0 , we expect to observe a mean of 200 molecules. The 
empirical CDF for each simulation model, compiled over 10^ 
realizations, is presented in Fig. The empirical CDFs of 
all three models match the Binomial CDF. Here, the Poisson 
approximation also matches the Binomial CDF, whereas the 































Table 2: Simulation test parameter summary. The 
units for k and Ativi are and ms, respectively. 


Test 

System 

Source 

Observer 

N 

k 

AIm 

1 

Fi 

Fi 

Half of Fi 

1000 

0 

5 

2 

All 

All 

Vi 

10800 

0 

10 

3 

All 

TX 

RX 

1500 

0 

5 

4 

- 

- 

- 

1000 

3 

Vary 

5 

All 

TX 

RX 

1500 

3 

5 



Figure 2: CDF of simulation 1, i.e., the uniform 
diffusion test where molecules are restricted to Vi 
only. The observation is made at time t = 5 s. 



4j= of Molecules Observed in Vi 


Figure 3: CDF of simulation 2, i.e., the uniform 
diffusion test that places molecules over all of V. The 
observation is made at time t = 20 s. 


Gaussian approximation is slightly less accurate. 

In the third simulation, we consider “point-to-point” diffu¬ 
sion, where 1500 molecules are released at the TX and then 
observed over time at the RX. The TX and RX are both 
squares of width 3 ym and are separated by a distance of 
15 ym from center to center, as shown in Fig. The time- 
varying channel impulse response, averaged over 10® real¬ 
izations, is plotted for different simulation models in Figs.j^ 
and[^ We omit a curve of the expected channel impulse re¬ 
sponse, as found by evaluating because it is effectively 
identical to that given by the MICRO simulation model. 

In Fig. 1^ we observe that the MESO model (with sub¬ 
volumes of size 3 ym everywhere) and the HYB model that 



Figure 4: Time-varying channel response of simula¬ 
tion 3, i.e., the “point-to-point” diffusion test. 


is microscopic in regions Vi and V 2 are both very accurate 
when compared with the MICRO model. The MESO-MS 
model with subvolumes of size 3 ym in V 2 is also very accu¬ 
rate, but increasing those subvolumes to 15 ym leads to an 
18 % overestimation of the channel impulse response at the 
time of the expected peak observation (~ 0.55 s). This ex¬ 
cess is because, from it takes longer for molecules to dif¬ 
fuse to larger subvolumes, which here are too close to the TX 
and RX. However, this model is still asymptotically accu¬ 
rate over time because the transition rate 0 derived in 
leads to a uniform molecule distribution. Finally, the HYB 
model that is only microscopic in Vi generally overestimates 
the channel impulse response (i.e., the MICRO/MESO in¬ 
terface is too close to the communication link). We further 
study the accuracy of HYB models in Fig. 

In Fig. we compare the MICRO model with variations 
of the HYB model that overestimated the channel impulse 
response in Fig. i.e., where regions V 2 and V 3 are both 
mesoscopic. Here, we vary the microscopic time step Atju 
of the HYB model, and we observe the resulting sensitiv¬ 
ity. Smaller time steps lead to fewer molecules at the RX 
(as expected; n\^2DAtM decreases but entering the MESO 
regime introduces the same uncertainty in a molecule’s lo¬ 
cation, leading to a net migration out of Vi). These results 
highlight the caution that must be taken when using hybrid 
models. We emphasize that we implemented simple transi¬ 
tion rules and that accuracy can be improved by optimizing 
the transition rules for a given time step as described in . 

In Fig.[^ we consider the empirical CDF of the third sim¬ 
ulation, evaluated at times t = 0.05 s and t — 0.2 s after 
the release by the TX, i.e., near the time of the peak of the 
expected signal and after the signal is expected to have de¬ 
creased by more than 3 dB from the peak value, respectively. 
At those times, from we expect 6.98 and 4.05 molecules, 
respectively. For clarity, we only plot the empirical CDFs 
for the MICRO model and the least accurate models pre¬ 
sented in Fig. since the CDFs for the MESO model with 
subvolumes of the same size and the HYB model with mi¬ 
croscopic V 2 are identical to that of the MICRO model. The 
MICRO model matches the Binomial CDF at both obser¬ 
vation times, and the Poisson approximation is effectively 
identical to the Binomial CDF. The simulation models that 
did not accurately capture the expected channel response 
also did not accurately match the Binomial CDF. Finally, 
we observe that the Gaussian approximation of the Bino- 













































Figure 5: Time-varying channel response of simula¬ 
tion 3, i.e., the “point-to-point” diffusion test, where 
we highlight the impact of Atm on the accuracy of 
the HYB model where /i 2 = /13 = 3 ^m. 



Figure 6: CDF of simulation 3, i.e., the “point-to- 
point” diffusion test. 


mial CDF is almost as poor as the least accurate simulation 
at both observation times, i.e., the MESO-MS model with 
/i 2 = 15 pm at t = 0.05 s and the HYB model that is meso¬ 
scopic in V 2 at t = 0.2 s. 

In the fourth simulation, we consider first-order degrada¬ 
tion only and do not allow molecules to diffuse. This test 
emphasizes the accuracy of simulating chemical reactions 
alone. We simulate the degradation of 1000 molecules when 
the reaction constant is A: = 3s“^. The time-varying re¬ 
sponse, averaged over lO'^ realizations, is observed in Fig. 
for the MICRO models with different values of Atm and the 
MESO model. We do not include a curve for the expected 
response, as given by ([^, but it is identical to the curve 
shown for the MESO model. We observe that the MICRO 
model is also very accurate for AAm varying over orders of 
magnitude, and the loss of accuracy when AIm = 50ms is 
only an artifact because this time step is longer than the 
observation period of 20 ms. 

In Fig. 1^ we observe the empirical CDF of the fourth 
simulation for the observation made t = 0.5 s after the start. 


Figure 7: Time-varying channel response of simula¬ 
tion 4, i.e., the first-order reaction test. Diffusion is 
not modeled. 



Figure 8: CDF of simulation 4, i.e., the first-order 
reaction test, where the observation is made at time 
t — 0.5 s. Diffusion is not modeled. 


From §, about 223 molecules are expected to remain. All 
of the simulation models agree with the Binomial CDF, and 
the Gaussian approximation is much more accurate than the 
Poisson approximation. Here, the underlying trial success 
probability is 0.223. 

In the fifth and final simulation, we combine the “point- 
to-point” diffusion test with first-order degradation. 1500 
molecules are released at the TX and then observed over 
time at the RX when the reaction constant is fe = 3s“^. 
Simulation results are averaged over 10® realizations. We 
observe that the accuracy of the simulations is consistent 
with that observed in the corresponding diffusion-only case. 
In Fig.[^ we observe the time-varying response for the same 
simulation models that we considered in the “point-to-point” 
diffusion test without degradation. We do not plot the ex¬ 
pected time-varying channel response, as given by the prod¬ 
uct of § and exp(—fct), because it is effectively the same 
as the average MICRO simulation. As in Fig. the MESO 
model, the MESO-MS model with /i 2 = 3 pm, and the HYB 
model where regions Vi and V 2 are microscopic yield aver¬ 
age simulation results that are very similar to the MICRO 
model. The remaining simulation models are noticeably less 
accurate than the MICRO model. 

In Fig. |10[ we observe the empirical CDF of the fifth sim¬ 
ulation for the observation made t = 0.05 s after molecules 
are released by TX. At that time, from and exp {—kt) or 
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Figure 9: Time-varying channel response of simu¬ 
lation 5, i.e., the ‘‘point-to-point” reaction-diffusion 
test. 
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Figure 10: CDF of simulation 5, i.e., the “point-to- 
point” reaction-diffusion test, where the observation 
is made at time t — 0.05 s. 


from inspection of Fig. about 6 molecules are expected. 
As in Fig.[^ we compare the least accurate simulation mod¬ 
els with the MICRO model, whose empirical CDF is once 
again equivalent to the Binomial CDF. Also, the Poisson 
approximation of the Binomial CDF is again much more 
accurate than the Gaussian approximation. 

5. CONCLUSIONS 

In this paper, we compared simulation models to assess 
their accuracy in diffusion, first order reaction, and first- 
order reaction-diffusion simulations. We observed the time- 
varying channel response and the empirical CDF at spe¬ 
cific time instants. The microscopic model and the meso¬ 
scopic model were generally accurate and their statistics 
very closely matched those of the underlying Binomial CDF 
for all simulations. We demonstrated that multi-scale and 
hybrid models could also maintain accuracy, unless we re¬ 
duced the computational complexity too close to the commu¬ 
nication link. Overall, the statistical accuracy of the receiver 
was not affected if the hybrid interface or transition to larger 
subvolumes was as far from the transmitter and receiver as 
the distance between the transmitter and receiver. 

We also compared the suitability of the Poisson and 
Gaussian approximations of the Binomial CDF, since these 
approximations are commonly applied in communications 
analysis. When a large fraction of the released molecules 


are expected at the receiver, the Gaussian approximation 
is more accurate. When a small fraction of molecules are 
expected, the Poisson approximation is more accurate. 

Our on-going work is the development of a molecular sim¬ 
ulator based on the models presented in this paper and the 
motivation in [^. Future implementation includes exten¬ 
sion to three dimensions, modeling fluid flow, and imple¬ 
menting more accurate rules for transitions between the mi¬ 
croscopic and mesoscopic regimes. 
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